Let ${\bf M}\in\mathbb{R}^{m_1\times m_2}$ be an unknown matrix with $r={\rm rank}({\bf M})\ll \min(m_1,m_2)$ whose thin singular value decomposition is denoted by ${\bf M}={\bf U}{\bf \Lambda}{\bf V}^{\top}$ where ${\bf \Lambda}={\rm diag}(\lambda_1,\cdots,\lambda_r)$ contains its non-increasing singular values. Low rank matrix regression refers to instances of estimating ${\bf M}$ from $n$ i.i.d. copies of random pair $\{({\bf X}, y)\}$ where ${\bf X}\in\mathbb{R}^{m_1\times m_2}$ is a random measurement matrix and $y\in\mathbb{R}$ is a noisy output satisfying $y={\rm tr}({\bf M}^{\top}{\bf X})+\xi$ with $\xi$ being stochastic error independent of ${\bf X}$. The goal of this paper is to construct efficient estimator (denoted by $\hat{\bf U}$ and $\hat{\bf V}$) and confidence interval of ${\bf U}$ and ${\bf V}$. In particular, we characterize the distribution of $$ {\rm dist}^2\big[(\hat{\bf U},\hat{\bf V}), ({\bf U},{\bf V})\big]=\|\hat{\bf U}\hat{\bf U}^{\top}-{\bf U}{\bf U}^{\top}\|_{\rm F}^2+\|\hat{\bf V}\hat{\bf V}^{\top}-{\bf V}{\bf V}^{\top}\|_{\rm F}^2. $$ We prove the asymptotical normality of properly centered and normalized ${\rm dist}^2\big[(\hat{\bf U},\hat{\bf V}), ({\bf U},{\bf V})\big]$ with data-dependent centering and normalization when $r^{5/2}(m_1+m_2)^{3/2}=o(n/\log n)$, based on which confidence interval of ${\bf U}$ and ${\bf V}$ is constructed achieving any pre-determined confidence level asymptotically.